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Abstract — For signals defined on the sphere, we develop 
a tool to project the signal onto the joint spatio-spectral 
domain. We revisit the definition of the spatially localized 
spherical harmonic transform (SLSHT) from the litera- 
ture, where azimuthally symmetric window functions were 
used for spatial localization. We propose a directional (i.e., 
azimuthally asymmetric) window function for spatial lo- 
calization, which enables the transform to reveal localized 
directional content in the spatio-spectral domain; hence, 
we call the proposed transform the directional SLSHT 
and the spatio-spectral representation as the directional 
SLSHT distribution. We present a fast algorithm for the 
efficient computation of the transform, which is required 
for practical applications. We study the accuracy and speed 
of our fast algorithm and show that it achieves very good 
numerical accuracy, with numerical errors at the level of 
fioating point precision. Finally, an example is presented 
to illustrate the capability of the directional SLSHT to 
resolve localized directional features in the spatio-spectral 
domain. 

Index Terms — Signal analysis, spherical harmonics, 2- 
sphere. 



L Introduction 

Signals that are inherently defined on the sphere 
appear in various fields of science and engineering, such 
as medical image analysis |[T|, geodesy ||2|, computer 
graphics ||3|, planetary science Q, electromagnetic in- 
verse problems f5l, cosmology [6], 3D beamforming fT\ 
and w^ireless channel modeling |8|. In order to analyze 
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and process signals on the sphere, many signal process- 
ing techniques have been extended from the Euclidean 
domain to the spherical domain (2), |[9|-|[23|. 

Due to the ability of wavelets to resolve localized 
signal content in both space and scale, wavelets have 
been extensively investigated for analyzing signals on the 
sphere ||9|, |[T3|-|[T5|, |[l9|-|[23j and have been utilized 
in various applications (e.g., in astrophysics p4|-[[29| 
and geophysics ^4J, pO| , PH). Some of the wavelet 
techniques on the sphere also incorporate directional 
phenomena in the spatial-scale decomposition of a sig- 
nal (e.g., [2T|-f23l). As an alternative to spatial-scale de- 
composition, spatio-spectral (spatial-spectral) techniques 
have also been developed and applied for localized spec- 
tral analysis, spectral estimation and spatially varying 
spectral filtering of signals |10l, pT2l, fTSl, fil\, f3^. 
The spectral domain is formed through the spherical 
harmonic transform which serves as a counterpart of the 
Fourier transform for signals on the sphere Q, |34[- 
|36). 

The localized spherical harmonic transform, composed 
of spatial windowing followed by spherical harmonic 
transform, was first devised in [TF] for localized spectral 
analysis. We note that the localized spherical harmonic 



transform was defined in |18| for azimuthally asym- 
metric (i.e., directional) window functions, however, it 
was applied and investigated for azimuthally symmet- 
ric functions only. Furthermore, a spectrally truncated 
azimuthally symmetric window function was used for 
spatial localization [ [T8| . Due to spectral truncation, the 
window used for spatial localization may not be concen- 
trated in the region of interest. This issue was resolved 
in p2| , where azimuthally symmetric eigenfunctions 
obtained from the Slepian concentration problem on 
the sphere were used as window functions (the Slepian 
concentration problem is studied for arbitrary regions on 
the sphere in f^). Following fTSf, the spatially localized 
spherical harmonic transform (SLSHT) for signals on 
the sphere has been devised in 1 10] to obtain the spatio- 
spectral representation of signals for azimuthally sym- 
metric window functions, where the effect of different 
window functions on the SLSHT distribution is studied. 
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Subsequently, the SLSHT has been used to perform 
spatially varying spectral filtering |12|, again with az- 
imuthally symmetric window functions. 

In obtaining the SLSHT distribution for spatio- spectral 
representation of a signal, the use of an azimuthally 
symmetric window function provides mathematical sim- 
plifications, however, such an approach cannot discrim- 
inate localized directional features in the spatio- spectral 
domain. This motivates the use of asymmetric window 
functions in the spatio-spectral transformation of a signal 
using the SLSHT. In order to serve this objective, we 
employ the definition of the localized spherical har- 



Concluding remarks are presented in Section |VII 



monic transform in [ [T8] | and define the SLSHT and 
the SLSHT distributions using azimuthally asymmetric 
window functions for spatial localization. Since the use 
of asymmetric window function enables the transform to 
reveal directional features in the spatio-spectral domain, 
we call the proposed transform the directional SLSHT. 
We also provide a harmonic analysis of the proposed 
transform and present an inversion relation to recover 
the signal from its directional SLSHT distribution. 

Since the directional SLSHT distribution of a signal 
is required to be computed for each spatial position and 
for each spectral component, and data- sets on the sphere 
are of considerable size (e.g., three million samples on 
the sphere for current data-sets |37| and fifty million 
samples for forthcoming data-sets [38]), the evaluation 
of the directional SLSHT distribution is computationally 
challenging. We develop fast algorithms for this purpose. 
Through experimental results we show the numerical 
accuracy and efficient computation of the proposed di- 
rectional SLSHT transform. Furthermore, due to the 
fact that the proposed directional SLSHT distribution 
depends on the window function used for spatial local- 
ization, we analyze the asymmetric band-limited window 
function with nominal concentration in an ellipsoidal 
region around the north pole, which is obtained from 
the Slepian concentration problem on the sphere. We 
also illustrate, through an example, the capability of 
the proposed directional SLSHT to reveal directional 
features in the spatio-spectral domain. 

The remainder of the paper is structured as follows. In 
Section [ll| we review mathematical preliminaries related 
to the signals on the sphere, which are required in the 
sequel. We present the formulation of the directional 
SLSHT, its harmonic analysis and signal reconstruction 
from the SLSHT distribution in Section [Till Different al- 
gorithms for the evaluation of the SLSHT distribution are 



provided in Section |IV| and the detailed analysis of the 
asymmetric window function is presented in Section |Vj 
In Section VI we show timing and accuracy results 
of our algorithms and an illustration of the transform. 



II. Mathematical Background 

In order to clarify the adopted notation, we review 
some mathematical background for signals defined on 
the sphere and the rotation group. 



A. Signals on the Sphere 

In this work, we consider the square integrable com- 
plex functions f{x) defined on unit sphere = {u G 
: |u| = 1}, where | • | denotes Euclidian norm, 
X = x{9,(j)) = (sin ^ cos 0, sin sin 0, cos 9)^ G 
is a unit vector and parameterizes a point on the unit 
sphere with 9 G [0,7r] denoting the co-latitude and 
(p G [0, 27t) denoting the longitude. The inner product 
of two functions / and h on is defined as [39J 



{f,h) 



f{x)h{x) ds{x) 



(1) 



where (•) denotes the complex conjugate, ds{x) = 
sin 9d9d(j) and the integration is carried out over the 
unit sphere. With the inner product in ([T]), the space of 
square integrable complex valued functions on the sphere 
forms a complete Hilbert space L^(S^). Also, the inner 
product in ([T]) induces a norm 11/11 A (/,/)i/2. We refer 
the functions with finite induced norm as signals on the 
sphere. 

The Hilbert space L^(S^) is separable and the spher- 
ical harmonics form the archetype complete orthonor- 
mal set of basis functions. The spherical harmonics, 
Yp{x) = Yp{9, (/)), for degree £ > and order |m| < £ 
are defined as ||5|, [ [36| 



YneA) ^NpPr (cos 9)e' 



where 



l2£ + l{£-my. 
Air {e + m)\ 



(2) 



(3) 



denotes the normalization constant and Pp are the 
associated Legendre polynomials defined, for m > 0, 
by: 



1) 



(4) 
(5) 



for \x\ < 1. With the above definitions, the spherical 
harmonics form an orthonormal set of basis functions, 
i.e., they satisfy {Yp^Yp') = See'Smm', where See' is 
the Kronecker delta. 
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By completeness and orthonormality of the spherical 
harmonics, we can expand any signal / G L^(S^) as 



(6) 



e=o m=-i 



where 



functions defined on the rotation group SO(3) and follow 
the orthogonality relation 

JS0(3) -h 1 

(11) 

where dp = dasin /3dl3dj and the integral is a triple 
integral over all rotations (a,/3,7) G SO(3) |36|. Thus, 
(/)^ A (^f^ yp) = / /(x)y^^(&)d5(&) (7) any function / G L^(SO(3)) may be expressed as 

/(/^) = E E E (12) 

where 

if)L,m' - ^ / /(P)^^^^'^P- (13) 

OTT JsO(3) 

The signal / is said to be band-limited with maximum 
degree Lf if {f)^^^, = 0, V£ > Lf. 



denotes the spherical harmonic coefficient of degree 
i and order m. The signal / is said to be band- 
limited with maximum spherical harmonic degree Lf if 
(/)7 = 0, V£>L/. 



B. Rotations on the Sphere and Wigner-D Functions 

Rotations on the sphere are often parameterized using 
Euler angles (q:,/3,7) G S0(3), where a G [0, 27r), 
/3 G [0, 7t] and 7 G [0, 27r) |36|. Using the 'zyz' 
Euler convention, we define the rotation operator 
for p = (a, /3, 7) G SO(3), which rotates a function on a 
sphere in the sequence of 7 rotation around z-axis, then 
/3 rotation about y-axis followed by a a rotation around 
z-axis. The spherical harmonic coefficient of a rotated 
signal ^pf is related to the coefficients of the original 
signal by 

i 

{%f)7 - E <.n'(p)(/)r''/^=(«'/3'7), (8) 

m'=-i 

where D^^^, (p) denotes the Wigner-L) function [ [36| of 
degree i and orders m and and is given by 

Di^^,ip) = Z)i^^,(a,/3,7) (9) 
= e-^-"<„,(/3)e-^-'^p 7), 

where d^^f{/3) is the Wigner-c? function |36| and is 
defined as 

n 

^(^ + my.{£ - my.{£ + m)\{£ - m)\ 
{l-\r m' — n)\{n)\{l — n — m)\{n — m' ^ m)\ 

2i—2n-\-m'—m / R\2n—m'-\-m 

sinf^] (10) 



COS 



2) 



C. Signals on the Rotation Group SO(3) 

For ^ > and m^m' ^ Z such that |m|, \m'\ < i, 
the Wigner-L) functions in ([9]) form a complete set 
of orthogonal functions for the space L^(SO(3)) of 



D. Discretization of and SO(3) 

In order to represent functions on and SO(3), it 
is necessary to adopt appropriate tessellation schemes 
to discretize both the unit sphere domain and the Eu- 
ler angle domain of SO(3). We consider tessellation 
schemes that support a sampling theorem for band- 
limited functions, which is equivalent to supporting an 
exact quadrature. 

For the unit sphere domain, we adopt the equiangular 
tessellation scheme |35| defined as (3^ = {One ^ 
^(2n^ + l)/(2L + l), cPn, = 2^V(2L + 1) : < < 
L, < < 2L}, which is a grid of (L + 1) x (2L + 1) 
sample points on the sphere (including repeated samples 
of the south pole) that keeps the independence between 
9 and for fix value of 9. For a band-limited function 
on the sphere / G L^(S^) with maximum spherical 
harmonic degree Lj, the sampling on the grid 
ensures that all information of the function is captured in 
the finite set of samples and, moreover, that exact quadra- 
ture can be performed [35 1. Note that this sampling 
theorem was developed only recently [^ | and requires 
approximately half as many samples on the sphere as 
required by alternative equiangular sampling theorems 
on the sphere [34 1. 

For the Euler angle representation of the rotation 
group SO(3), we consider the equiangular tessellation 
scheme (Bl 2Tina/{2L + 1), 

2Tin^/{2L + 1), 7^^ = 27rn^/(2L + 1) : < n^,,/!^ < 
2L, < < L}. Again for a function / G L2(S0(3)) 
with maximum spectral degree Lf, the sampling of a 
function / on €7,/ ensures that all information of the 
function is captured and also permits exact quadrature 
(which follows from the results developed on the sphere 



|35|). 
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III. Directional SLSHT 

We describe in this section the directional SLSHT, 
which is capable of revealing directional features of 
signals in the spatio-spectral domain. For spatial localiza- 
tion, we consider the band-limited azimuthally asymmet- 
ric window function which is spatially concentrated in 
some asymmetric region around the north pole. Since the 
rotation around the z-axis does not have any affect on an 
azimuthally symmetric function, the localized spherical 
harmonic transform using an azimuthally symmetric 
window function can be parameterized on the sphere. 
However, if an azimuthally asymmetric window is used 
to obtain localization in the spatial domain, the rotation 
of the window function is fully parameterized with the 
consideration of all three Euler angles (a, ^5, 7) G SO(3). 
We refer to the spatially localized transform using an 
asymmetric window as the directional SLSHT. Here, 
we first define the directional SLSHT distribution which 
presents the signal in the spatio-spectral domain. Later in 
this section, we present the harmonic analysis of SLSHT 
distribution and provide an inversion relation to obtain 
the signal from its given directional SLSHT distribution. 

A. Forward Directional SLSHT 

Definition 1 (Directional SLSHT): For a signal / G 
L^(S^), define the directional SLSHT distribution com- 
ponent g{p;i, m) G L^(SO(3)) of degree i and order m 
as the spherical harmonic transform of a localized signal 
where localization is provided by the rotation operator 
acting on window function h G L^(S^), i.e.. 



f{x){%h){x)Yp{x)ds{x) (14) 



for < 



< Lg, \m\ 



< i, where Lg Lf -\- L^ 
denotes the maximum spherical harmonic degree for 
which the distribution components g{p;i^m) are non- 
zero, and Lf and L^ denote the band-limits of the 
signal / and the window function h, respectively. Also, 
each distribution component g{p;i,m) is band-limited 
in p = (a,/3,7) G S0(3) with maximum degree L^, 
i.e., when expressed in terms of Wigner-D functions. We 
elaborate this shortly. Furthermore, we consider unit en- 
ergy normalized window functions such that (/i, h) = 1. 
Remark 1: The directional SLSHT distribution com- 



ponent in ( [T4| ) can be interpreted as the spherical 
harmonic transform of the localized signal where the 
window function h provides asymmetric localization at 
spatial position = a) G and the first rotation 
determines the orientation of the window function at 
X. If the window function is azimuthally symmetric, 
this orientation of the window function by 7 becomes 



invariant and the SLSHT distribution components are 
defined on L2(§2) ||io|. 

Since the maximum spectral degree for which the 
SLSHT distribution is defined is = + L/^, we 
consider the window function such that L^ < Lf to 
avoid extending Lg significantly above Lf. We discuss 
the localization of the window function in spatial and 
spectral domains later in the paper. 



B. Harmonic Analysis 

We now present the formulation of the directional 
SLSHT distribution if the signal / and the window 
function h are represented in the spectral domain. Using 
the expression of the spherical harmonics of a rotated 
function in ([8]), we can write the SLSHT distribution 
component g{p;i,m) in ( [T4| ) as 



p=0 q=—p q' = —p 

where 

T{l',m'-p,q-l,m)^ [ Yir\x)Y^^{x)Yl^ds{x) 

(16) 

denotes the spherical harmonic triple product, which 
can be evaluated using Wigner-3j symbols or Clebsch- 
Gordan coefficients [36], |40|. 

Remark 2: By comparing g{p; m) in ( fTS) ) with ([12]), 
we note that the band-limit of g{p;i^m) in p is given 
by Lh. Since < Lf and p < L^ in ( [T5] ), our statement 
that the distribution component g{p; m) is non-zero for 
^ ^ Lg = Lf + Lh follows since the triple product in 



( [T6| ) is non-zero for £ < Lf + Lh only. 



C. Inverse Directional SLSHT 

Here, we define the inverse directional SLSHT to 
reconstruct a signal from its SLSHT distribution. The 
original signal can be reconstructed from its direc- 
tional SLSHT distribution through the spectral domain 
marginal, that is, by integrating the SLSHT distribution 
components over the spatial domain SO(3) |18|. Using 
our harmonic formulation in ( [T5] ), define (/)^^^ as the 
integral of the SLSHT distribution component g{p; m) 
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over SO(3) giving 



(/): 



SO(3) 



g{p;i,m)dp, 0<i<Lf 



E E (/)"EE E W^' 

i'=Om' = -i' p=Oq=-pq' = -p 

xTie',m';p,q;£,m) [ Dl^,{p)dp 



SO(3) 



(17) 



where we have used the following property of Wigner-Z) 
functions 



/ 

JSO(3) 



DlAp)dp = 8^2^5,0 5,^0, (18) 



which stems from ( [TT] ). Using the expression in ( fTT) ), we 
can find the spherical harmonic coefficient (/)^ of the 
signal / as 



if)] 



(19) 



which indicates that we only need to know the DC 
component of the window function (/i)^ in order to 
obtain the signal from its directional SLSHT distribution. 
It further imposes the condition that the DC component 
of the window function must be non-zero. Although the 
distribution components in ( [T5] ) are defined up to degree 
Lg = Lf -\- Lh, we only require the components up to 
Lf for signal reconstruction. 

Since the directional SLSHT distribution components 
g{p;i^m) in ([14]) are defined for 



< Lg, the number 
2 while 



of distribution components are of the order Lg 
the sampling of p is of the order thus, the di- 
rect evaluation of the directional SLSHT distribution 
is prohibitively computationally expensive. Therefore 
efficient algorithms need to be developed which reduce 
the computational complexity. We address this problem 
in the next section. 

IV. Efficient Computation of Directional 
SLSHT Distribution 

Here, we present efficient algorithms for the com- 
putation of the directional SLSHT distribution of a 
signal and the signal reconstruction from its directional 
SLSHT distribution. First, we discuss the computational 
complexities if the SLSHT distribution components are 
computed using direct quadrature as given in ([14]) or us- 
ing the harmonic formulation in ([15]). Later, we develop 
an alternative harmonic formulation which reduces the 
computational burden. Finally, we present an efficient 



algorithm that incorporates a factoring of rotations [ |4T| 
and exploits the FFT. 

First we need to parameterize the required tessellation 
schemes for for the representation of the signal / and 
the window h and for SO(3) which forms the spatial 
domain of the directional SLSHT distribution. Since the 
maximum spectral degree of the signal f is Lf, we 
therefore consider the equiangular tessellation to 
represent signals on L^(S^). Since the maximum degree 
for all SLSHT distribution components g{p; m) in p is 
Lh, we therefore consider the tessellation (Bl^ to repre- 
sent the SLSHT distribution components on i7^(SO(3)). 



A. Direct Quadrature and Harmonic Formulation 

We define the forward spatio-spectral transform as 
evaluation of each SLSHT distribution component 
g{p;i,m). Evaluation of the forward spatio-spectral 



transform using exact quadrature in ([14]) requires the 
computation of two dimensional summation over the 
tessellation of for each 3-tuple (a,/3,7). Since there 
are 0{L^ such 3-tuples in the tessellation scheme 
and the SLSHT distribution components are of the 
order 0{L'j), the computational complexity to compute 
all distribution components using direct quadrature is 
0{Lj^Ll). Using the harmonic formulation in ( fTS) ), the 
complexity to compute each SLSHT distribution compo- 
nent is 0{L'jLfJ and to compute all SLSHT distribution 
components is O(LjL^). Although the harmonic formu- 
lation in ( [T5] ) is useful to establish that the signal can be 
reconstructed from the directional SLSHT distribution, 
it is much more computationally demanding than direct 
quadrature. We develop efficient algorithms in the next 
subsection which improve the computational complexity 
of the harmonic formulation and make it more efficient 
than direct quadrature. 

For the inverse directional SLSHT distribution, we 
only need to integrate over SO(3) to obtain the signal 
in the spherical harmonic domain as proposed in ([17]). 
Since the integral can be evaluated by a summation over 
all Euler angles using quadrature weights, an efficient 
way to recover the signal from its SLSHT distribution 
is through direct quadrature, with complexity of 0{LfJ 
for each distribution component and 0{L'jLl) for all 
components. 

1) Exact Quadrature and Quadrature Weights: In 
order to evaluate the integral in ( [TT] ) exactly, we need 
to define quadrature weights along Euler angle /3 in the 
tessellation (Bl^- We evaluate the integral in ( [TT] ) by the 
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2Lu 



2L, 



following summation 

^^^^ ^{2Lh + lf ^n^o^o 

tIq, =0 72/3 =U =L) 

'yn^;£,m)q{/3ni,), (20) 

where the quadrature weights qiPm,) are evaluated by 
employing the method outlined in (351, giving 



4.r2(L^J + l)-\ 



/3n3 =0 



Stt^ ^ w{—m) cosm/3^^, otherwise 

(21) 



where it;(m) is defined as |35| 

m = ±l, 

w{m) = ^ 0, modd, m 7^ 1, (22) 



l-m2 ' 



meven. 



= E E E (/^D; 



Lh P P 

E 

7?=0 q=—p q' = —p 



Using the spherical harmonic transform method in 



|35|, [hy can be computed in 0{LfJ. The product 



inside the inner summation of ([24]) can be computed in 
0{LI) for each p and for all < p. The inner 



two summations of ([24]) include complex exponentials 
and therefore we can employ an FFT to evaluate these 
two summations in 0(Lf^\og2 L^) for each p and each 
/3. Since, both p and /3 are of the order L^, each 
SLSHT distribution component can be computed using 
([24]) in 0(L^ min(log2L/,L/,) + L^logsL/,) and all 



components can be computed with overall complex- 
ity of 0(L^ min(log2L/,L/,) + L^L^ logs L/,). It is 
important to note that we have considered the factor 
min(log2 L/, L^), which arises from the computation 
of the spherical harmonic coefficients of fY^, in the 
total complexity and this may dominate the other factor 
L'jLf^\og2 Lh if Lh <^ Lf. In any case, the alterna- 
tive harmonic formulation in ([24]) is more efficient to 



B. Forward Directional SLSHT: Semi-Fast Algorithm 

We consider an alternative harmonic formulation in 
order to reduce the computational complexity of the 
forward SLSHT. We may write the directional SLSHT 
distribution component g{p;i^m) in ([14]) as a spherical 
convolution [T5 1 of h and the spherical harmonic mod- 
ulated signal fYp, giving 

9{p-e,m) = E E E WWX{hy^Dl^'i^,P,^), 

p=0 q=—p q' = —p 

(23) 

which can be expressed, using the definition of the 
Wigner-D function in ([9]), as 



compute than the harmonic formulation in ( [T5] ) or direct 
quadrature of ([14]). 

Using the factoring of rotations approach which was 
first proposed in |41| and has been used in the imple- 
mentations of the fast spherical convolution [ |42J and 
the directional spherical wavelet transform p3|, the 
complexity of the harmonic formulation in ( [24] ) can be 
reduced further. We present this development in the next 
subsection. 



C. Forward Directional SLSHT: Fast Algorithm 

Here we employ a factoring of rotations to develop 
an algorithm which reduces the factor L'jLf^ log2 Lh in 
the complexity of the semi-fast case and thus improves 
the performance of the overall algorithm. It also reduces 
the asymptotic complexity of the overall algorithm when 
0{L'iLf^\og2 Lh) dominates in the overall computa- 



tional complexity of the formulation in (24), i.e., when 



X <,.(/3)e-^^'^e-^^". (24) ^ logs Lh > Lj min(log2 Lf, Lh). 



The band-limit of the spherical harmonic modulated 
signal / is Lf-\-£. Since the maximum £ for which g 
is non-zero is Lf+Lh, we must compute up to / Y^^^j^^, 
which is band-limited to 2Lj + L/^. However, we only 
need to compute the spherical harmonic coefficients 
(^fYf^y^ of the modulated signal up to degree p < Lh- 
Therefore, the computation of spherical harmonic trans- 
form of fYp is an interesting sub-problem. By either 
using the separation of variable approach or the method 
outlined in [35j , the spherical harmonic coefficients 
(7^r)p forO <p< Lh, \q\ < p of the signal jYl^ 
can be computed in 0{L'j min(log2 L/, Lh)) time for a 
given i and m (see Appendix[A] for details). 



By factoring the single rotation by (a,/3,7) into two 
rotations (IS), |[4l|. 



Pl = {a- 7T/2, -7T/2, P2 = (0, 7T/2, 7 + 7T/2), 

(25) 

and noting the effect of rotation on spherical harmonic 
coefficients in ([8]), we can write the Wigner-L) function 
in ([9]) as 

q" = -p 

(26) 



7 



where A^^, = ^/ (^/2) and we have used the following 
symmetry properties of Wigner-rf functions [40] 

= = (27) 

Using the Wigner-D expansion given in ([26]), we can 
write the alternative harmonic formulation of the SLSHT 
distribution component g{p;i^m) in ( [23] ) as 

p=0 q=—p q' = —p 

X ^l'q^l','e-"i''-"i"^-'^'\ (28) 

q" = -p 

where p = (a,;5,7). By reordering the summations we 
can write 

Lh Lh Lh 

g{p;£,m) ^ J2 Yl Yl '^q,Q',Q"(^^^) 

q=-Lh q'=-Lh q"=-Lh 

^^-iqa-iq"P-q'T^ p= (0,^,7), (29) 

where 

p=max(|g|,|g'|,|g"|) 

(30) 

Comparatively, the computation of the SLSHT distribu- 
tion components using the expression given by ( [29] ) is not 
more efficient than the initial expression ( [24] ). However, 
the presence of complex exponentials can be exploited 
by employing FFTs to evaluate the involved summations. 

1) Efficient Implementation: The objective of factor- 
ing the rotations is to carry out the /3 rotation along the 
?/-axis as a rotation along the z-axis. The rotations along 
the z-axis are parameterized using complex exponentials 
and thus these rotations can be applied with much less 
computational burden, by exploiting the power of an 
FFT, relative to a rotation about the ?/-axis. All the three 
rotations which characterize the spatial domain of the 
SLSHT distribution components appear in complex ex- 
ponentials in ( [29] ) and thus we can use FFTs to evaluate 
the summation of Cq^q' ^q"{l^m) over q, q' and q" . First 
we need to compute Cq^q' ^q"[l^m) for each I and for 
each m which requires the one-dimensional summation 
over three dimensional grid formed by q, q' and q" and 
thus can be computed in 0{Lf^). Using Cq^q' ^q"[l^m), 
the summation over the complex exponentials in ( [29] ) 
can be carried out in 0{L^\og2Lh) using FFTs. The 
overall complexity of this approach is dominated by the 
computation of Cq^q' ^q"[l^m), that is, 0{Lj^ for each 
SLSHT distribution component and 0{L'iLfJ for the 



complete SLSHT distribution. We note that the evalu- 
ation of Cq^q' ^q"{l^m) rcquircs the computation of A^^, 
which can be evaluated over the (g, q') plane for each 
p using the recursion formula of [43] with a complexity 
of 0{Lf^). Since p is of the order L^, the A matrices 
can be evaluated in 0{LfJ, which does not change 
the overall complexity of our proposed algorithm. The 
overall asymptotic complexity of our fast algorithm is 
thus 0{Ljmm{log2Lf,Lh) + L}Ll). 

Remark 3: If we want to analyze the signal / 
with multiple window functions, then we do not need 
to recalculate the spherical harmonic transform of 
the modulated signal fY^, which accounts for the 
0(Lj min(log2 L/, L^)) factor in the overall complexity. 
Once it is computed, the SLSHT distribution can be 
computed in O(LjL^) time for each window function 
of the same band-limit using the proposed efficient 
implementation. 

Our proposed formulation and efficient implementa- 
tion can be further optimized in the case of a steerable 
window function. Steerable functions have an azimuthal 
harmonic band-limit in m that is less than the band-limit 
in i (see [20] for further details about steerability on the 
sphere). In this case, the L'jLf^ factor contributing to the 
overall asymptotic complexity of the fast algorithm is 
reduced to L'jLf^- Furthermore, we may then compute 
the directional SLSHT for any continuous 7 G [0, 27r) 
from a small number of basis orientations (due to the 
linearity of the SLSHT). 

If the signal and window function are real, the compu- 
tational time can be further reduced by considering the 
conjugate symmetry relation of the spherical harmonic 
coefficients. Furthermore, in this setting, the SLSHT 
distribution components also satisfy the conjugate sym- 
metry property 

g{p;£, -m) = (-1)^ ^(p; ^, m) (31) 

and we do not need to compute the SLSHT distribution 
components of negative orders. 

V. Window localization in Spatial and 
Spectral Domains 

The directional SLSHT distribution is the spherical 
harmonic transform of the product of two functions, the 
signal / and the rotated window function h and we 
must be careful in interpreting the directional SLSHT 
distribution in the sense that we do not mistake using 
the signal to study the window because there is no 
distinction mathematically. The window function should 
be chosen such that it provides spatial localization in 
some spatial region around the north pole (origin). Since 
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we have considered a band-limited window function, 
the window function cannot be perfectly localized in 
the spatial domain due to the uncertainty principle on 
the sphere [44J . However, it can be optimally localized 
by maximizing the energy concentration of the window 
function in the desired directional region |2|. 

The interpretation and the effectiveness of the di- 
rectional SLSHT distribution depends on the chosen 
window function. The window function with maximum 
localization in some defined asymmetric region provides 
directional localization and thus reveals directional fea- 
tures in the spatio- spectral domain. The more directional 
the window function, the more directional features it can 
reveal in the spatio- spectral domain but the larger the 
maximum spherical harmonic degree L^. Recall that the 
maximum degree of the directional SLSHT distribution 
components is given by Lg = Lf -\- L^. Thus, when 
the signal is expressed in the spatio- spectral domain its 
spectral domain is extended by L^, which results in 
spectral leakage. Therefore, we want the window func- 
tion to be simultaneously maximally localized in some 
spatial region ^ C and have the minimum possible 
band-limit which achieves the desired level of energy 
concentration in the spatial region This problem is 
well known in time-frequency analysis [ |45| and has been 
extended to spatio- spectral analysis on the sphere for 
azimuthally symmetric window functions [ [T0| , where the 
minimum band-limit azimuthally symmetric eigenfunc- 
tion obtained from the Slepian concentration problem 
on the sphere p2| has been proposed as a suitable 
window function fTO|. 

Here, we propose using a band-limited eigenfunction 
obtained from the solution of the Slepian concentration 
problem [2| as a window function, concentrated in a 
spatially localized ellipsoidal region around the north 
pole. We first parameterize the ellipsoidal region and 
later analyze the resulting eigenfunctions from the per- 
spective of the uncertainty principle on the sphere |44|. 
We note that the choice of the asymmetric region as an 
ellipsoidal region is only one possibility and the analysis 
can be extended to other asymmetric regions such as strip 
regions around the north pole [ ,46 J . 

A. Parametrization of Window Function 

We consider a band-limited window function of max- 
imum spherical harmonic degree L/^, which is spatially 
concentrated in the asymmetric ellipsoidal region on the 
sphere with major axis along the x-axis, and thus is 
orientated along the x-axis. The ellipsoidal region can 
be parameterized using the focus 9^ of the ellipse along 
the positive x-axis and the arc length a of the semi-major 



axis: 



+ A,((0,(/)),(0e,^)) <2a}, 



(32) 



where < 6>c < 7r/2 and 9^ < a < 7r/2. Here 
A5((6>,(/)), (6>',(/)0) = arccos (sin6>sin(9'cos((/) - 00 + 
cos cos 0^ denotes the angular distance between two 
points (0,0) and {9' ^(j)') on the sphere. 

Remark 4: For a given focus 9^, the region becomes 
more directional as the arc length a approaches 9c. For 
a = 7r/2, the region becomes azimuthally symmetric, 
i.e., we recover the polar cap of central angle 7r/2. 
As an illustration, the ellipsoidal regions ^[o^^a] with 
focus 6>c = 7r/6 and a G {tt/G + 7r/24o"7r/6 + 
7r/120, 7r/6 + tt/SO} are shown in Fig.[l] 

B. Slepian Concentration Problem 

As a result of the Slepian concentration problem |[2|, 
|46| to find the band-limited function with bandwidth 
Lh and maximal spatial concentration in an ellipsoidal 
region ^(51^^), we obtain {L^ + 1)^ eigenfunctions. The 
eigenvalue associated with each eigenfunction serves as 
a measure of the energy concentration in the spatial 
region. Here we consider the use of the band-limited 
eigenfunction with maximum energy concentration in the 
ellipsoidal region for given band-limit and refer to 
such an eigenfunction as the eigenfunction window. If A 
denotes the area of the ellipsoidal region, it is shown ||2|, 
[ [46| that most of the eigenvalues lie either near zero or 
unity for both symmetric and asymmetric regions and the 
sum of all eigenvalues, referred as an equivalent of the 
Shannon number |2], is equal to A^o = ^(^/i + l)^/(47r). 
Also, it is shown empirically in [ [32| that there exist less 
than or equal to A'o — 1 spatially concentrated eigen- 
functions with non-insignificant energy concentration. 
By noting these developments and empirical results, and 
considering that must be chosen that we obtain at 
least one eigenfunction which is spatially concentrated in 
the ellipsoidal region, we recover the following empirical 
lower bound for L^'. 




A 



(33) 



where [ (•) ] denotes the integer ceiling function (this 
bound follows directly from A'o ^ 2). We analyze this 
bound later in this section. 

Let Amax denote the eigenvalue associated with the 
most spatially concentrated eigenfunction. By finding 
the minimum value of Lh which ensures that Amax is 



(a) (b) 
Fig. 1: Ellipsoidal regions ^(o^^a) on the sphere with focus 9c 

a = 7r/6 + Tv/UO and (c) a = 7r/6 + 7r/240. 
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= 7t/6 and major axis: (a) a = tt/G + tt/SO, (b) 
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(a) (b) (c) 

Fig. 2: Band-limited eigenfunction windows h obtained from Slepian concentration problem on the sphere with 
maximal spatial concentration in an ellipsoidal region ^(51,, a) of focus 9c = 7r/6 and major axis: (a) a = 7r/6+7r/80, 
(b) a = 7r/6 + 7r/120 and (c) a = tt/G + 7r/240. 



greater than or equal to the desired energy concentra- 
tion, an eigenfunction window h{x) with desired energy 
concentration in the ellipsoidal region and minimum 
possible band-limit can be obtained. Thus, the focus 
of an ellipsoidal region 9c, arc length of semi-major axis 
a and maximum spectral degree fully parameterize 
the eigenfunction window. As an illustration, the three 
eigenfunction windows having 90% spatial concentration 
in the ellipsoidal regions ^(o^^a) with focus 9c = tt/Q 
and a G {tt/G + 7r/240, tt/G + ^120, tt/G + tt/SO} and 
respective maximum spherical harmonic degree G 
{18, 14, 11} are shown in Fig.|2] For these ellipsoidal 
regions, in Fig.|3] we show A^ax versus spherical har- 
monic band-limit L denoting the band-limit of the most 
concentrated eigenfunction window, where we note the 
difference in and the bound for the desired 
concentration level. For 90% desired concentration level, 
we observe that the spherical harmonic degree which 
ensures Amax > 0.9 deviates further from the bound 



Bh given by ( [33] ) as the concentration region becomes 
more directional (as expected since the bound does not 
incorporate the level of directionality of the region). 

C. Analysis of Eigenfunction Window Concentrated in 
Ellipsoidal Region 

1) Directionality: For the use of the eigenfunction 
window in obtaining the directional SLSHT distribution, 
the directionality of the eigenfunction window must be 
a key criterion in selecting the eigenfunction window 
for the identification of localized features of a signal in 
the spatio-spectral domain. Also, the directionality of the 
eigenfunction window depends on the directionality of 
the ellipsoidal region. We use the definition of the auto- 
correlation function on the sphere as a measure of the 
directionality of the eigenfunction window |20|. Since 
the major axis and minor axis of an ellipsoidal region lie 
along X-axis and y-axis respectively, the auto-correlation 
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Fig. 3: The maximum eigenvalue A^ax of the eigen- 
function h with spherical harmonic band-limit L for 
spatial concentration in an ellipsoidal region of focus 
9c ^ "k/Q and major axis a as indicated. Black vertical 
lines indicate the bound on the band-limit given 



by ( [33] ) and the gray vertical lines indicate the actual 
Lh which ensures 90% concentration in an ellipsoidal 
region. 



function is defined as the inner product of eigenfunction 
window with its version rotated around the z-axis: 

Ch{<;) = {Dp,h, /i) , PI = (0,0,0, ^ ^ [0, |], (34) 

which can be expressed in the harmonic domain as 

Lh e 



^im(; \j^m\2 



(35) 



£=0 m= 



Due to the symmetry of the ellipsoidal region and the 
eigenfunction window, we have considered in the range 
of to 7r/2. 

For the eigenfunction windows presented in the pre- 
vious subsection, the auto-correlation function for each 
window is shown in Fig.|4| which indicates that Ch{^) 
decays more rapidly from unity at c;^ = for more 
directional window and therefore, the peakedness of 
Chi^) quantifies the ability of the window function to 
reveal the localized directional features of a signal in 
the spatio-spectral domain. 

2) Spatial and Spectral Localization: Here we study 
the spatial and spectral localization of eigenfunction win- 
dows from the perspective of the uncertainty principle on 
the sphere [ [44| , [ [47| , according to which, the function 
cannot be simultaneously localized in both the spatial 
and spectral domains. The following inequality specifies 
the uncertainty principle for unit energy functions de- 
fined on the sphere, which relates the trade-off between 
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Fig. 4: Auto-correlation function Ch{^) as described in 



( |34| ) for eigenfunction windows h with spatial concen- 
tration in an ellipsoidal region of focus 6>c = and 
major axis a as indicated. 
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Fig. 5: Spatial variance (200-3), spectral variance {aO 



and uncertainty product, as described in d36b, (37) and 



( [38] ), are shown for eigenf unctions with spatial concen- 
tration in an ellipsoidal region of focus 9^ ^ n/Q and 
major axis a as indicated. We have plotted a scaled value 
of spatial variance (20cr|) for better visualization. 



the spatial and spectral localization of a function: 

• (JL > 1, 



(36) 



1 



at 



where as and ctl denote the variance of the band-limited 
unit energy eigenfunction window in the spatial domain 



and spectral domain respectively and are defined as |47| 



1 - 




sm 



{29) \h{9,^)\^ d9d^ 



(37) 
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and 



m\2 
t I • 



(38) 



^=0 



Due to the consideration of unit energy functions, < 
(JS < I. We note that smaller variance indicates better 
localization of the window function. The variance in 
the spatial and spectral domains and the uncertainty 
product for eigenfunctions concentrated in an ellipsoidal 
region of focus tt/G and major- axis of different values 
are shown in Fig. [5] As expected, the variance in the 
spatial domain decreases as the region becomes more 
directional, whereas the variance in the spectral domain 
increases because increases. We also note that the 
uncertainty product increases as the region becomes 
more directional. 

VI. Results 

In this section, we first demonstrate the numerical 
validation and computation time of our algorithms to 
evaluate the directional SLSHT components. Later, we 
provide an example to illustrate the capability of the 
directional SLSHT, showing that it reveals the directional 
features of signals in the spatio- spectral domain. The 
implementation of the our algorithms is carried out in 
MAT LAB, using the MAT LAB interface of the SSHt[^ 
package (the core algorithms of which are written in 
C and which also uses the FFTV\|^ package to compute 
Fourier transforms) to efficiently compute forward and 



inverse spherical harmonic transforms |35|. 



A. Numerical Validation and Computation Time 

In order to evaluate the numerical accuracy and the 
computation time, we carry out the following numerical 
experiment. We use the band-limited function h for 
spatial localization with band-limit L/^ = 18 and spatial 
localization in the region ^(7r/6,7r/6+7r/240)- We generate 
band-limited test signals with band-limits 18 < < 
120 by generating spherical harmonic coefficients with 
real and imaginary parts uniformly distributed in the 
interval [0, 1]. 

For the given test signal, we measure the computation 
time Ti to compute all directional SLSHT distribution 
components g{p;i^m) for i < Lf + Lh and m < \i\ 
using our fast algorithm presented in Section |IV-C[ 
where we compute the Wigner-d functions on-the-fly 
for the argument 7r/2 by using the recursion of Tra- 
pani [|43|. We also record the time T2 to compute the 



spherical harmonic transform of the modulated signal 
only and the time T3 to recover a signal from its SLSHT 
distribution components. For the computation of the 
spherical harmonic transform of the modulated signal, 
we use the separation of variable approach presented in 
Appendix |Aj All numerical experiments are performed 
using MAT LAB running on a 2.4 GHz Intel Xeon pro- 
cessor with 64 GB of RAM and the results are averaged 
over ten test signals. The computation time ti and T2 
in seconds are plotted against the band-limit Lf of the 



test signal in Fig. 



6a 



which scale as 0{Lj^) and thus 



corroborates the theoretical complexity. We also note that 
the major cost in the evaluation of SLSHT distribution 
components is the computation of the spherical harmonic 
transform of the spherical harmonic modulated signal. 
The computation time ts_ for the inverse directional 
SLSHT is plotted in Fig 



6b 



which evolves as 0{L^ 



for 



fixed Lh, again corroborating the theoretical complexity. 

We reconstruct the original signal from its SLSHT 
distribution components using ([20]) and ( fTT) ), in order 
to assess the numerical accuracy of our algorithms by 
measuring the maximum absolute error between the 
original spherical harmonic coefficients of the test signal 
and the reconstructed values. The maximum absolute 



error is plotted in Fig. 6c for different band-limits Lf, 
which illustrates that our algorithms achieve very good 
numerical accuracy with numerical errors at the level of 
floating point precision. 

B. Directional SLSHT Illustration 

In order to illustrate the capability of the proposed 
transform to reveal the localized contribution of spec- 
tral contents and probe the directional features in the 
spatio-spectral domain, we analyze the Mars topographic 
map (height above geoid) as a signal on the sphere, 
which is obtained by using the spherical harmonic model 
of the topography of Mar^ The Mars topographic map 
is shown in Fig. [7] in the spatial domain, where the grand 
canyon Valles Marineris and the mountainous regions 
of Tharsis Montes and Olympus Montes are shown, 
leading to the high frequency contents. We note that 
the mountainous regions are non-directional features of 
the Mars map, whereas the grand canyon serves as a 
directional feature with direction orientated along a line 
of approximate constant latitude. The spherical harmonic 
coefficients certainly provide details about the presence 
of higher degree spherical harmonics in the signal, 
but do not reveal any information about the localized 
contribution of higher degree spherical harmonics. 



^ http://www.jasonmcewen.org/ 




^ http://www.fftw.org/ 




' http://www.ipgp.fr/~wieczor/SH/ 
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(a) Forward transform computation time (b) Inverse transform computation time (c) Numerical validation 

Fig. 6: Numerical validation and computation time of the proposed algorithms, (a) The computation time ti (in 
seconds) is the time to compute all SLSHT distribution components of the test signal with band-limit Lf and is 
plotted for test signals of different band-limits. The computation time T2 denotes the time to compute the spherical 
harmonic transform of the modulated signal. Both ti and T2 evolve as 0(Lj) as shown by the solid red line, (b) 
The computation time T3 denotes the time to compute the inverse directional SLSHT, which scales as 0{L'j) (for 
fix Lh) as shown by the solid red line, (c) The maximum error e between the spherical harmonic coefficients of the 
test signal and those of the signal obtained via reconstruction from the SLSHT distribution components is plotted, 
which empirically appears to scale as 0{L), as shown by the solid red line. 




Fig. 7: Mars signal in the spatial domain. The grand 
canyon Valles Marineris and the mountainous regions of 
Tharsis Montes and Olympus Montes are indicated. 



If we analyze the signal by employing the SLSHT 
using an azimuthally symmetric window function, the 
presence of localized contributions of higher degree 
spectral contents can be determined in the spatio- spectral 
domain p^, however the presence of directional features 
cannot be extracted. Here, we illustrate that the use 



of the directional SLSHT enables the identification of 
directional features in the spatio- spectral domain, which 
is due to the consideration of an asymmetric window 
function for spatial localization. 

We obtain the directional SLSHT distribution com- 
ponents g{p;i^m) of the Mars map / using the band- 
limited eigenfunction window h with = 60 and 90% 
concentration in the spatial domain in an ellipsoidal re- 
gion ^(7r/i2,77r/i2)- The SLSHT distribution components 
g{p;i^m) for order m = 15 and degrees 20 < ^ < 27 
and 50 < ^ < 57 are shown in Fig. [8] and Fig.|9] respec- 
tively, where the components in the panels are for Euler 
angle (a) 7 = and (b) 7 = 1007r/201 = L5630 ^ 7r/2. 
Since the ellipsoidal region is oriented along the x-axis, 
the window with orientation 7 = provides localization 
along colatitude and the window with orientation 7 ^ 
7t/2 provides localization along longitude. It is evident 
that using orientation of the window 7 ^ 7r/2 probes the 
information about the grand canyon (directional feature) 
along longitude in the spatio- spectral domain. The local- 
ized contribution of higher degree spherical harmonics 
towards the mountainous region can also be observed 
in Fig.[8] for degree 20 < ^ < 27 and both 7 = and 
7 = 7r/2. However, there is not significant contribution 
of spherical harmonics of degree 50 < ^ < 57 towards 
mountainous region as indicated in Fig. [9^ where 7 = 0, 
but the localization of the directional features along the 
orientation 7 ^ 7r/2 are revealed in the spatio- spectral 
domain as shown in Fig.|9|3. Due to the ability of the 
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Fig. 8: Magnitude of the components of the directional SLSHT distribution of the Mars signal obtained using the 
eigenfunction window concentrated in an ellipsoidal region of focus 9c = 7r/12 and major axis a = In/ 12. For 
fixed orientation 7, the distribution components g(p;i^m) are mapped on the sphere using p = (0,6^,7) for order 
m = 15 and degree 20 < i < 27. The components are shown for orientation (a) 7 = and (b) 7 ^ 7r/2 of the 
window function around the z-axis. Top left: g{p; 20, 15), top right: g{p; 23, 15). 



directional SLSHT to reveal the localized contribution of 
spectral contents and the directional or oriented features 
in the spatio-spectral domain, it can be useful in many 
applications where the signal on the sphere is localized 
in position and orientation. 

VII. Conclusions 

We have presented the directional SLSHT to project a 
signal on the sphere onto its joint spatio-spectral domain 
as a directional SLSHT distribution. In spirit, the direc- 
tional SLSHT is composed of SO(3) spatial localization 



followed by the spherical harmonic transform. Here, we 
have proposed the use of an azimuthally asymmetric 
window function to obtain spatial localization, which 
enables the transform to resolve directional features 
in the spatio-spectral domain. We have also presented 
an inversion relation to synthesize the original signal 
from its directional SLSHT distribution. Since data- 
sets on the sphere are of considerable size, we have 
developed a fast algorithm for the efficient computation 
of the directional SLSHT distribution of a signal. The 
computational complexity of computing the directional 
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Fig. 9: Magnitude of the components of the directional SLSHT distribution of the Mars signal obtained using the 
eigenfunction window concentrated in an ellipsoidal region of focus 9c = 7r/12 and major axis a = In/ 12. For 
fixed orientation 7, the distribution components g(p;i^m) are mapped on the sphere using p = (0,6^,7) for order 
m = 15 and degree 50 < ^ < 57. The components are shown for orientation (a) 7 = and (b) 7 ^ 7r/2 of the 
window function around z-axis. Top left: g(p; 50, 15), top right: g{p; 53, 15). 





SLSHT is reduced by providing an alternative harmonic 
formulation of the transform and then exploiting the 
factoring of rotation approach |41 1 and the fast Fourier 
transform. The computational complexity of the pro- 
posed fast algorithm to evaluate SLSHT distribution of 
a signal with band- limit L f using window function with 
band-limit is 0(Lj min(log2 L/, L/j) + L^'f^t) 
compared to the complexity of direct evaluation, which 
is 0{Lj^LfJ. The numerical accuracy and the speed of 
our fast algorithm has also been studied. The directional 
SLSHT distribution relies on a window function for spa- 



tial localization; we have analyzed the band-limited win- 
dow function obtained from the Slepian concentration 
problem on the sphere, with nominal concentration in 
an ellipsoidal region around the north pole. We provided 
an illustration which highlighted the capability of the 
directional SLSHT to reveal directional features in the 
spatio- spectral domain. 

There are natural extensions and applications of the 
work presented. Our directional SLSHT distribution can 
be used in many applications to resolve localized di- 
rectional contents in a signal. Furthermore, the use of 
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the directional SLSHT transform to perform spatially- 
varying filtering of the directional features of a signal 
in the spatio- spectral domain needs to be studied |T2l. 
Since the SLSHT is analogous in nature to the short- 
time Fourier transform (STFT) in time-frequency analy- 
sis f48l, following the extension of the STFT presented 
in [491 , JSOl , the use of different window functions 
for obtaining localization at different spatial positions 
presents an open problem for future work. 

Appendix A 
Spherical Harmonic Transform of 
Modulated Signal 

Our objective is to compute the spherical harmonic 
transform of the modulated signal fYf^, for given i and 
m, up to degree L^, which can be performed in two 
ways. Firstly, we could use a simple separation variable 
technique given by 



27T 



fie,' 



sm0d0. (39) 



Since < i < Lf + L^, we need to consider the 
modulated signal fY^ sampled on the grid ©2L/+2L^ 
for the explicit evaluation of exact quadrature f35l. Using 
([39]), the integral over (p can be computed first in 0{L'j) 
for each q and in 0{L'jLh) for all \q\ < L^. We can 
also employ an FFT to evaluate the integral over (f) in 
0(Ljlog2L/) time. Then, using the exact quadrature 
weights derived in |[35|, the integral over 9 can be 
computed in 0{Lf) for each q and for in 0{LfLf^) 
for all q. Thus the overall complexity using this 
approach is 0(max(min(Lj log2 L/, L'jL^)^ LfLfJ) ^ 
0(Lj min(log2 L/, L/z)), since throughout we assume 
Lh < Lf. Alternatively, we can employ a factoring 
of rotation approach to compute the spherical har- 
monic transform according to the method given in [ 35J , 
which has the complexity 0(max(Lj log2 L/, LfJ) ^ 
0{L'j log2 Lf) but supports a sampling of the modulated 
signal on the lower resolution grid &2Lf+Lh' Therefore, 
the spherical harmonic transform of the modulated signal 
fYl^ up to degree can be computed with asymptotic 
complexity 0(L'j min(log2 L/, Lh)). 
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